#!/usr/bin/python
'''
Comptues validation statistics for an assembly.
'''
import sys
import fasta_read
import numpy as np

# parameters.
sam_file = sys.argv[1]
valid_fasta = sys.argv[2]

# build reference.
print "reading reference."
ref_data = fasta_read.read_fasta_dict(valid_fasta)

tmp = {}
for name in ref_data:
	n = name.split(" ")[0]
	tmp[n] = ref_data[name]
ref_data = tmp

# build reference.
print "building reference."
ref_cov = {}
for name in ref_data:
	# make numpy.
	ref_cov[name] = np.zeros(len(ref_data[name]), dtype=np.int) 

# read sam file.
print "reading sam file."
fin = open(sam_file, "rb")
lines = fin.readlines()
fin.close()

# parse coverage.
print "parsing sam file."
for line in lines:
	# tokenize.
	if line[0] == "@": continue
	tmp = line.split("\t")
	start = int(tmp[3])
	ref = tmp[2].replace(">","")
	cigar = tmp[5]
	
	# check we have ref.
	if ref not in ref_cov:
		print "bad ref"
		print ref
		sys.exit()
		
	# parse cigar.
	clist = []
	num = ""
	for c in cigar:
		# determine if number.
		try:
			i = int(c)
		except:
			i = None
	
		# not number.
		if i == None:
			# add pair.
			clist.append((c, int(num)))
			num = ""
			continue
			
		# is a number.
		num += c
			
	# Add hits to reference.
	sz = len(ref_cov[ref])
	i = start
	for pair in clist:
		# skip non-matches.
		if pair[0] != "M":
			i += pair[1]
			continue
		
		# note match.
		if pair[1] >= sz:
			stop = sz
		else:
			stop = pair[1]
		for j in range(i,stop):
			ref_cov[ref][j] += 1
		i += pair[1]
		
# parse out percent aligned.
print "results"
bins = [.99, .95, .80, .50, 0]
cnts = [0, 0, 0, 0, 0]
z = len(bins)
for name in ref_cov:
	# Get total.
	tot_hit = np.sum(ref_cov[name] > 0)
	tot = ref_cov[name].size
	pl = float(tot_hit) / float(tot)

	for i in range(z):
		if pl > bins[i]:
			cnts[i] += 1
	
	#print name, tot_hit, tot
	
# pull out statistics.
for j in range(z):
	print "%f\t%i" % (bins[j], cnts[j])

	

